Visualizing the data stream of long-term records for the online book

Load the R packages

library(ggplot2) # graphing
library(tidyr) # dataframe manipulation
library(lubridate) # for handling date times
library(plotly) # for making interactive graphs

Read in data from EDI repository

# Package ID: knb-lter-hbr.59.10 Cataloging System:https://pasta.edirepository.org.
# Data set title: Hubbard Brook Experimental Forest: Daily Temperature Record, 1955 - present.
inUrl1  <- "https://pasta.lternet.edu/package/data/eml/knb-lter-hbr/59/10/9723086870f14b48409869f6c06d6aa8" 
infile1 <- tempfile()
try(download.file(inUrl1,infile1,method="curl"))
if (is.na(file.size(infile1))) download.file(inUrl1,infile1,method="auto")


air <-read.csv(infile1,header=F ,skip=1 ,sep="," , 
               col.names=c("date",     
                 "STA","MAX","MIN","AVE","Flag"),check.names=TRUE)
unlink(infile1)

# We have a data column, but its not formatted as a date
air$DATE<-ymd(air$date) # change how R interprets Date to be a date
air$Year<-year(air$DATE)

packageID<-"Package ID: knb-lter-hbr.59.10"

Calculate annual air temp for each water year

# Use water year starting June 1st
w.year <- as.numeric(format(air$DATE, "%Y"))
june.july.sept <- as.numeric(format(air$DATE, "%m")) < 6
w.year[june.july.sept] <- w.year[june.july.sept] - 1
air$wyear<-w.year # add water year as a column to precip dataset

#subset to just HQ
hq<-air[air$STA=="HQ",]

# make sure you are only using complete years for the record
pre.obs<-as.data.frame(table(hq$STA , hq$wyear))
pre.obs[pre.obs$Freq<360, "Use"]<-"incomplete wyear" # incomplete is < 360 days
pre.obs[is.na(pre.obs$Use),"Use"]<-"complete"

# match by water year-station, only use 'complete' years
pre.obs$wys<-paste(pre.obs$Var2, pre.obs$Var1)
hq$wys<-paste(hq$wyear, hq$STA)
hq$Use<-pre.obs$Use[match(hq$wys, pre.obs$wys)]
hq.complete<-hq[hq$Use=="complete",]

# aggregate to annual average
hqa<-aggregate(list(MAX=hq.complete$MAX, MIN=hq.complete$MIN, AVE=hq.complete$AVE), by=list(wyear=hq.complete$wyear), FUN="mean")

hq_g<-gather(hqa, "type","value", 2:4)

## calculate difference between recent 6 years compared to first 6 years
magmin<-format(round(mean(tail(hqa$MIN))-mean(head(hqa$MIN)), digits=0))
magave<-format(round(mean(tail(hqa$AVE))-mean(head(hqa$AVE)), digits=0))
magmax<-format(round(mean(tail(hqa$MAX))-mean(head(hqa$MAX)), digits=0))

# order temp types
hq_g$type<-factor(hq_g$type, levels=c("MAX","AVE","MIN"))

graph air temperature at the Hubbard Brook Headquarters

## make a ggplot graph
a1<-ggplot(hq_g, aes(x=wyear, y=value, col=type))+geom_point()+ geom_line()+
  ylab("Air temperature (C)")+xlab("Water year (June 1st)")+
  theme_bw()+theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
  labs(col='')+
  scale_color_manual(values=c(MIN="blue",AVE="black",MAX="orange"))+
  geom_smooth(method="lm", se=F)+
  theme(axis.text.y=element_text(size=15, colour = "black"),
        axis.title.y=element_text(size=17, colour = "black"),
        axis.text.x=element_text(size=15, colour = "black"))+
   geom_text(aes(x=1975, y=14.5, label = paste0("Average daily high ",magmax, " degree hotter")),
        size = 4.5, color="orange")+
   geom_text(aes(x=1972, y=8.2, label = paste0("On average ",magave, " degrees hotter")),
        size = 4.5, color='black')+
   geom_text(aes(x=1975, y=3, label = paste0("Average daily low ",magmin, " degrees hotter")),
        size = 4.5, color='blue')+
  ggtitle("Air temperature at the Hubbard Brook Experimental Forest Headquarters")

## plotly object
p1<-ggplotly(a1)
                                  
p1

Write a html file

# this line writes the html file to create interactive graphs for the online book
htmlwidgets::saveWidget(as_widget(p1), "HQ_air_temp.html")

hubbardbrook.github.io/HQ_air_temp.html